#include "stdafx.h"
#include "ckSIMProcessor.h"
#include "ckRasterBand.h"

namespace ck
{
	namespace algorithm
	{
		namespace raster
		{
			SIMProcessor::SIMProcessor()
			{
				m_pBand1 = NULL;
				m_pBand2 = NULL;
				m_sim = 0.0;
				m_k1 = 0.01;
				m_k2 = 0.03;
			}


			SIMProcessor::~SIMProcessor()
			{
				if (m_pBand1)
				{
					m_pBand1->Release();
					m_pBand1 = NULL;
				}

				if (m_pBand2)
				{
					m_pBand2->Release();
					m_pBand2 = NULL;
				}
			}

			RESULTCODE SIMProcessor::Execute()
			{
				if (!m_pBand1 || !m_pBand2)
				{
					return RC_FAILURE;
				}

				z_uint32 w1, h1, w2, h2;
				w1 = m_pBand1->GetWidth();
				h1 = m_pBand1->GetHeight();
				w2 = m_pBand2->GetWidth();
				h2 = m_pBand2->GetHeight();
				if ((w1 != w2) || (h1 != h2))
				{
					return RC_FAILURE;
				}

				ck::data::ckPixelType type1, type2;
				type1 = m_pBand1->GetPixelType();
				type2 = m_pBand2->GetPixelType();
				if (type1 != type2)
				{
					return RC_FAILURE;
				}

				z_uint32 pixelSize = m_pBand1->GetPixelSize();
				double L = pow(2, (pixelSize)) - 1;

				double c1 = pow(m_k1 * L, 2);
				double c2 = pow(m_k2 * L, 2);

				ck::data::ck_raster_stats stats1, stats2;
				m_pBand1->ComputeStatistics(&stats1);
				m_pBand2->ComputeStatistics(&stats2);

				double m1 = stats1.mean;
				double m2 = stats2.mean;

				double var1 = pow(stats1.stdev, 2);
				double var2 = pow(stats2.stdev, 2);
				double coVar = CalculateCoVariance(m_pBand1, stats1.mean, m_pBand2, stats2.mean);

				m_sim = (2 * m1*m2 + c1) * (2 * coVar + c2) / ((m1*m1 + m2 * m2 + c1)*(var1 + var2 + c2));

				

				return RC_SUCCESS;
			}

			void SIMProcessor::SetBand(data::RasterBand* pBand1, data::RasterBand* pBand2)
			{
				m_pBand1 = pBand1;
				m_pBand2 = pBand2;
				if (m_pBand1)
				{
					m_pBand1->AddRef();
				}
				if (m_pBand2)
				{
					m_pBand2->AddRef();
				}
			}

			double SIMProcessor::GetValue()
			{
				return m_sim;
			}

			double SIMProcessor::CalculateCoVariance(data::RasterBand* pBand1, double mean1, data::RasterBand* pBand2, double mean2)
			{
				double coVar = 0.0;

				ck::data::ckPixelType type = pBand1->GetPixelType();
				switch (type)
				{
				case ck::data::ckPixelTypeByte:
					coVar = CalculateCoVariance_Byte(pBand1, mean1, pBand2, mean2);
					break;
				case ck::data::ckPixelTypeUInt16:
					coVar = CalculateCoVariance_UInt16(pBand1, mean1, pBand2, mean2);
					break;
				}

				return coVar;
			}

			double SIMProcessor::CalculateCoVariance_Byte(data::RasterBand* pBand1, double mean1, data::RasterBand* pBand2, double mean2)
			{	
				double sum = 0;
				z_uint32 count = 0;
				z_uint32 w = pBand1->GetWidth();
				z_uint32 h = pBand2->GetHeight();

				pBand1->ReadData();
				pBand2->ReadData();
				
				z_byte* ptr1 = (z_byte*)pBand1->GetData();
				z_byte* ptr2 = (z_byte*)pBand2->GetData();				

				z_uint32 i, j;
				for (i = 0; i < h; i++)
				{
					for (j = 0; j < w; j++, ptr1++, ptr2++)
					{
						//if ((*ptr1) == noData || (*ptr2) == noData)
						//{
						//	continue;
						//}

						sum += (*ptr1 - mean1) * (*ptr2 - mean2);
						count++;
					}
				}

				double coVar = ((double)sum) / (count - 1);

				pBand1->ReleaseData();
				pBand2->ReleaseData();

				return coVar;
			}

			double SIMProcessor::CalculateCoVariance_UInt16(data::RasterBand* pBand1, double mean1, data::RasterBand* pBand2, double mean2)
			{
				double sum = 0;
				z_uint32 count = 0;
				z_uint32 w = pBand1->GetWidth();
				z_uint32 h = pBand2->GetHeight();

				pBand1->ReadData();
				pBand2->ReadData();

				z_uint16* ptr1 = (z_uint16*)pBand1->GetData();
				z_uint16* ptr2 = (z_uint16*)pBand2->GetData();

				z_uint32 i, j;
				for (i = 0; i < h; i++)
				{
					for (j = 0; j < w; j++, ptr1++, ptr2++)
					{
						//if ((*ptr1) == noData || (*ptr2) == noData)
						//{
						//	continue;
						//}

						sum += (*ptr1 - mean1) * (*ptr2 - mean2);
						count++;
					}
				}

				double coVar = ((double)sum) / (count - 1);

				pBand1->ReleaseData();
				pBand2->ReleaseData();

				return coVar;
			}
		}
	}
}


